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Abstract 

Using data from e + e~ annihilation into hadrons, taken with the OPAL detec- 
tor at LEP at the Z pole between 1991 and 1995, we performed a simultaneous 
measurement of the colour factors of the underlying gauge group of the strong in- 
teraction, Cf and Ca, and the strong coupling, a s . The measurement was carried 
out by fitting next-to-leading order perturbative predictions to measured angular 
correlations of 4-jet events together with multi-jet related variables. Our results, 

C A = 3.02 ± 0.25(stat.) ± 0.49 (syst.) , C F = 1.34 ± 0.13(stat.) ± 0.22(syst.) , 

a s (M z ) = 0.120 ± 0.011(stat.) ± 0.020(syst.) , 

provide a test of perturbative QCD in which the only assumptions are non-abelian 
gauge symmetry and standard hadronization models. The measurements are in 
agreement with SU(3) expectations for Cf and Ca and the world average of a s {Mz)- 
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1 Introduction 



Electron-positron annihilation into hadrons at high energies provides a precise means 
to test Quantum Chromodynamics (QCD), the theory of strong interactions. Due to 
the purely leptonic initial state, there are many experimental quantities for which the 
long-distance (non-perturbative) effects are modest, especially when compared to similar 
quantities in hadron-hadron collisions or deep inelastic scattering. These quantities, for 
instance the total cross section and jet-related correlations, can be calculated in perturba- 
tive QCD as a function of a single parameter, the strong coupling strength, a s . Therefore, 
many QCD tests based on measurements of a s have been carried out at LEP, the Large 
Electron- Positron collider at CERN (see, for example, [0-0 )■ 

A key ingredient of QCD is the underlying gauge group. The simultaneous measure- 
ment of the strong coupling and the eigenvalues of the quadratic Casimir operator of 
the underlying gauge theory, the Cp and Ca colour factors, provides a more general and 
comprehensive test of QCD than measurements of a s alone. The possible existence of a 
light gluincQ influences both a s and the measured value of the colour factors (or, assuming 
SU(3) dynamics, the number of light fermionic degrees of freedom Nf). A simultaneous 
fit of these parameters to data provides a means to investigate whether the data favour 
or exclude the additional degrees of freedom. 

Previous analyses of LEP data have been performed to probe the underlying gauge 
structure of the theory |9|-[T3|. The leading order perturbative prediction was fitted to 
distributions of four-jet angular correlations (normalised to the total number of the se- 
lected four-jet events), leaving the ratios of the colour factors as free parameters. Due 
to the normalization, the leading order prediction of these distributions is independent of 
the strong coupling. However, the independence of the leading order prediction from the 
strong coupling does not imply that the radiative corrections are negligible. In fact, in 
Ref. jnj the value of T R /Cp (see the Appendix) was estimated to increase about 25% if 
next-to-leading order (NLO) theoretical predictions are used instead of the leading order 
(LO) ones. It is therefore desirable to explicitly verify the effect of the next-to-leading 
order corrections on these normalised angular distributions. 

Beyond leading order in the theoretical predictions, the strong coupling and the colour 
factors are interdependent. On the one hand, even the normalised next-to-leading order 



predictions for the angular distributions depend upon a s p^j , while on the other hand, 
a s itself depends on the colour factors through its running. In this sense, beyond leading 
order, the only consistent way to determine the three free parameters is a simultaneous 
measurement of the strong coupling and the colour factors. 

Predictions beyond leading order are made possible by theoretical developments achieved 
in the last few years. For multi-jet rates as well as numerous event shape distributions 
with perturbative expansions starting at 0(a s ), matched next-to- leading order and next- 
to-leading logarithmic approximations provide very precise descriptions of the data over 
a wide range of the available kinematic region P^oL]!^]. Also, recently, the next-to-leading 
order predictions for the distributions of four-jet angular correlations have been calcu- 



lated |TJ|, rS|. To make use of these developments, we perform a simultaneous fit of the 
strong coupling and ratios of colour factors using next-to-leading order predictions of 



1 The exclusion of gluinos with mass between 3-5 GeV has been debated in the literature recently 
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four-jet rates, differential two-jet rates and four-jet angular correlations. The data were 
collected by the OPAL Collaboration at LEP. 

The outline of the paper is as follows: In Section 2, we present the observables used 
in the analysis and describe their best available perturbative predictions. The analysis 
procedure is explained in detail and the fit results are presented in Section 3. Section 4 
contains the discussion of the systematic checks which were performed and the corre- 
sponding systematic errors. We collect our results and compare them to the results of 
previous studies in Section 5. Section 6 contains our conclusions. 



2 Observables 

To perform a simultaneous measurement of the strong coupling and the colour factors we 
use multi-jet related variables to gain sensitivity to a s and four-jet angular correlations 
to gain sensitivity to the colour factors. 

There are many ways of defining jets. For the present analysis we use the Durham 



scheme |16|,|20l. Starting by defining each particle to be an individual jet, a resolution 



variable y^ is calculated for each pair of jets i and j: 

2 mia(Ef, E]) (1 -cos0, 



Vij 



e 



VIS 



where Ei and Ej are the energies of jets i and j, 6ij is the angle between them and E V - 1S is 
the sum of the energies of the visible particles in the event or of the partons in a theoretical 
calculation. If the smallest value of is smaller then a predefined value y C ut, the pair 
is replaced by a pseudo-jet with four momentum p^ = Pi + Pj and the clustering starts 
again with momenta and pj dropped and p^ added to the final state. Clustering ends 
when the smallest value of y^ is larger than y cut . 

In our analysis we use the differential two-jet rate, 

„ / 1 da , , 

D 2 2/23 = -r— , 2 

O"tot d?/23 

where 7/23 is the y cut value for which the two- and three-jet configurations are separated 
in a given event, and the four-jet rate, 

R,(y C ut) = ff4 ' jct(?/cut) , (3) 

Otot 

to gain sensitivity to a s . For theory, <r tot is the total hadronic cross section at 0(a s ) 
accuracy. For experiment, cr to t is the total visible hadronic cross section. The D2 dis- 
tribution has been extensively used at LEP to determine a s [|,2l|. Use of R4 allows us 



to obtain sensitivity to a s using a four-jet based variable. For these observables both 
next-to-leading order (NLO) P~7|, |T5|] and next-to-leading logarithmic (NLL) [[LlJ pertur- 



bative predictions are known. The former provide a fixed order approximation valid for 
large values of the resolution variable, while the latter provide an all order summed but 
approximate cross section valid for small values of the resolution variable. The explicit 
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NLO and NLL formulae and their matched expressions used in our fits are given in the 
Appendix. Using those formulae, we obtain various theoretical predictions for the four-jet 
rates (R4) and the differential two-jet rates (D 2 ) which are shown in Fig. [1|. In these plots, 
the theoretical predictions with the exception of the fitted curves were obtained using the 
world average of the strong coupling a s (Mz) = 0.118 |22[ and standard QCD values for 
the colour factors, Cp = 4/3 and Ca = 3. The theoretical predictions in Fig. [I] are shown 
in comparison to OPAL data corrected to the parton level. The manner in which the data 



are corrected to the parton level is explained in Section |3.3| . It is seen that the theoretical 
predictions provide a good description of the data. For the theoretical plots in Fig. [I], 
we used the i?-matched expression for i? 4 and lni?-matched expression for D 2 (see the 
Appendix). 

The other class of observables we employ in our measurements is the four-jet angular 
correlations. The angular variables are 



the Bengtsson-Zerwas angle xbz = ^[(i>i x P2), {.Pz x Pa)]', 



the modified Nachtmann-Reiter angle f24|: Qnr = ^[(Pi ~ P2), (P3 — Pi)] 



the Korner-Schierholtz-Willrodt angle 
®ksw = \ UiiPi x P4), (p 2 x p 3 )] + Z[(pi x p 3 ), (p 2 x p 4 )]) ; 



the angle between the two lowest energy jets flOf : a 3i 4 = l[p3,p±], 



where the denote the three- momenta of the energy-ordered jets (E 1 > E 2 > E 3 > E^). 
We defined the jets using the Durham clustering and selected four-jet events at y cut = 
0.008. In order to have a high statistics four-jet sample, y cut should be chosen in the 
range from 10 -2 to 10~ 3 . However, we cannot choose y cut to be much smaller than 0.01 for 
technical reasons (see Section 3.3.1). For the four-jet events, used in obtaining the angular 
correlations, an additional cut was placed on the energy of each jet: E > E niin = 3 GeV. 

To obtain as much information as possible, we used all four angular variables in our 
analysis. The shape of the distributions of these angular correlations at LO and NLO are 
hardly distinguishable. Fig. shows the normalised NLO theoretical predictions fitted to 
the data. 

Although the various angular correlations are not entirely independent, the correla- 
tions (see Section O ) are not large as can be seen from Table |I|, where the maximum 



values of the magnitudes of the bin-wise correlations between the observables are given. 
We also performed fits using a single angular variable together with one of the jet-related 
variables and we observed that these fits give very scattered central values. However, in 
these cases we found very strong correlations between the fitted parameters (~ 1), indi- 
cating that the distributions do not contain sufficient information to constrain all of the 
QCD parameters simultaneously. 



3 Analysis Procedure and Results 
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COS 0NR cos $ksw cos a 34 


cos xbz 

COS ©NR. 
COS $KSW 


0.46 0.34 0.45 
0.37 0.37 
0.44 



Table 1: Maximum values of the bin-wise correlations between the angular correlation 
variables. 

3.1 The OPAL detector 



A detailed description of the OPAL detector can be found in Ref . |26| . This analysis relies 
mainly on the reconstruction of charged particle trajectories and on the measurement of 
energy deposited in the electromagnetic and hadronic calorimeters. 

The central detector contains a silicon micro-vertex detector and three drift chamber 
devices: an inner vertex chamber, a large jet chamber and a surrounding ^-chamber. The 
central detector is located inside a solenoidal magnet which provides a uniform axialQ 
magnetic field of approximately 0.435 T. Tracking of the charged particles is performed 
with this detector. Most of the tracking information is obtained from the jet chamber 
which provides up to 159 measured space points per track, and close to 100 % tracking 
efficiency in the region |cos#| < 0.98. The average angular resolution is approximately 
0.1 mrad in <p and 1 mrad in 6. 

Electromagnetic energy is measured by lead glass calorimeters surrounding the magnet 
coil, separated into a barrel (|cos0| < 0.82) and two end-cap (0.81 < | cos 6*| < 0.98) 
sections. The electromagnetic calorimeter consists of 11 704 lead glass blocks with a 
depth of 24.6 radiation lengths in the barrel and more than 22 radiation lengths in the 
end-caps. This calorimeter is surrounded by the hadronic calorimeter of the sampling 
type measuring the energy of hadrons emerging from the electromagnetic calorimeter. 

3.2 Data samples 

The analysis presented in this paper is based on 4.1 million hadronic events recorded 
within 3 GeV of the Z peak by the OPAL detector between 1991 and 1995. The descrip- 
tions of the OPAL trigger system and the offline multihadronic event selection are given 



in Refs. |27| and ||28|| . In this analysis, additional criteria were applied to eliminate poorly 
measured tracks and obtain well contained events. 

Each track was required to have (i) at least 40 measured points in the jet chamber, (ii) 
transverse momentum in the r—<f> plane greater than 0.15 GeV/c, (iii) measured momentum 
less than 60GeV/c, (iv) | cos 6 1 ! < 0.97, (v) distance of closest approach to the origin 
in the r—<f> plane of no more than 5 cm, (vi) and less than 25 cm in the z direction. 
Hadronic events were required to contain at least five tracks to reduce contamination 
from e + e~ — > t + t~ events. 



2 In the OPAL coordinate system the +x axis points towards the centre of the LEP ring, the y axis 
points approximately upwards and the z axis points in the direction of the electron beam. The polar 
angle 9 and the azimuthal angle <fi are defined with respect to z and x, respectively, while r is the distance 
from the z-axis. 
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Clusters in the electromagnetic calorimeter were accepted if they had more than 
0.1 GeV energy in the barrel section or more than 0.25 GeV in the end-cap section. The 
corresponding clusters were required to span at least two lead glass blocks. Clusters in 
the hadronic calorimeter were required to have energies larger than 1 GeV. 

For the calculation of the visible energy, E v i S , the tracks are assumed to have the pion 
mass, and clusters are treated as photons. To determine the energy of each cluster mea- 
sured in the electromagnetic and hadronic calorimeters, we used a matching algorithm |2{| . 
This algorithm reduces double counting of energy and gives a better resolution both in 
energy and angle than the calorimeters alone. We did not impose energy-momentum 
conservation, but checked that such a constraint does not affect our results. 

The event thrust axis was determined using the accepted tracks and clusters. In order 
that the event be well contained, we required | cos thrust! < 0.9. 

The final event sample contained about 3.6 million hadronic Z decays from which 
about 250 000 were identified as four-jet events. 



3.3 Data Correction 

We choose as the theoretical reference the distributions given in Eqs. (0), ( p8| ) of the 
Appendix which are obtained from a parton level calculation. In order to compare our 
(detector level) data to these (parton level) theoretical predictions the measured distri- 
butions were corrected for the effects of the detector and hadronization. 

First, we combined the distributions of the six variables cos xbz, cos Onr, cos $ksw ; 
cos«34, D2 and R4 into a one-dimensional distribution, with 128 bins, 20 for each angular 
correlation and 24 each for the four-jet rates (covering the range of 0.001 < y cu t < 1.0) 
and the differential two-jet rates (in the range of 1.0 < — ln?/23 < 5.8). Next, we applied 
bin-by-bin corrections for detector distortions and hadronization, as described below. The 
correction factors were obtained from the same distributions computed from Monte Carlo 
events. 



3.3.1 Monte Carlo simulation 

The correction factors were estimated using two different modes of the JETSET Monte 
Carlo program |3^,|3l]|: the parton shower mode and the matrix element mode. In the 
parton shower approach, the evolution of the parton system is treated as a branching 
process based on the leading logarithmic approximation starting from the original qq 
pair. In the matrix element mode, up to four partons are generated according to the 
second-order ERT calculation p2[ . 

It is well known that the parton shower versions of the available Monte Carlo programs 
provide a good description of jet rates, but cannot describe the structure of the four-jet 
events properly ||33|| . Therefore, to determine correction factors for the four-jet angular 



correlations, we used Monte Carlo events with only four partons (qqgg and qqqq) in the 
final state. These events were generated using the second-order matrix element (ME) 
calculation, subjected to hadronization, as implemented in JETSET version 7.4 assuming 
standard QCD colour factors and the set of parameters ERT-MC-1 of Ref. [S3]. 560 000 



Monte Carlo events were generated. To avoid singularities in the four-parton generation 



9 



of the LO matrix elements, an intrinsic cutoff yo has to be applied in the Monte Carlo 
simulation. In JETSET, y = 0.01 is used, calculated according to the invariant mass res- 
olution definition, which for four-jet events corresponds to a Durham resolution of about 
0.004 (yielding the same four-jet rates). Our resolution criterion of y cut = 0.008 is safely 
larger than the intrinsic cut-off. For normalised angular correlations, the distributions at 
leading and next-to-leading order are almost the same [14,19]. Thus for the purpose of 
determining the corrections the ERT Monte Carlo events are sufficient. 

For the correction of the differential two-jet rates and the four-jet rate, we used 2 mil- 
lion events generated with the parton shower (PS) version of JETSET, using parameters 
tuned to OPAL data at y/s =91 GeV described in Ref. []35|| . For both models, hadroniza- 



tion of the parton system is based on the Lund string- fragmentation scheme |36|]. The 



events of both types of samples were passed through the OPAL detector simulation |3j 
processed by the same reconstruction programs, and subjected to the same event selection 
criteria as the data. 



3.3.2 Correction procedure 

To correct the data for the finite acceptance and resolution of the detector, we determined 
bin-by-bin correction factors for every bin of the distribution from Monte Carlo samples 
with hadronization included, both before (hadron level, Hi MC ) and after (detector level, 
.Dj MC ) the detector simulation and selection cuts. The hadron level treats all charged 
and neutral particles with lifetimes greater than 3 x 10 -10 s as stable and has no initial 
state photon radiation. The correction factor for detector effects for the i th bin of the 
distribution is given by: 

TT MC 

r det - — (a\ 

We examined the effects of initial state photon radiation by switching the radiation 
on in the hadron level Monte Carlo simulation. The correction factors were found to be 
identical within statistical uncertainties with or without photon radiation over the entire 
phase space. 

We estimated the hadronization effects by computing the distributions of Monte Carlo 
events at the parton (Pj MC ) and hadron levels (Hi MC ) (before and after hadronization) 
and defined the correction factor for hadronization by: 

p.MC 

(i had _ £j /r\ 
* ~~ rj MC • \°) 



Using the two correction factors Cj et and Cj a , we correct the measured distribution 
to the parton level according to 

Di COTT = d det C7> ad A meas • (6) 

We show the distributions of the total correction factors, Cj tot = d det Ci had , in Figs. (TJ 
and |[ 
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3.4 Fit procedure 



Having prepared the corrected distribution of six variables, we performed a \ 2 minimiza- 
tion to determine the most probable values of the variables rj = 0£ s Cf/{2tt), x = Ca/Cf 
and y = Tr/Cf (see Section [A]) with the program MINUIT [38|. The function \ 2 was 
defined by: 



X 2 = E« 6i 8. 



where 8i is the bin-wise difference between the theoretical QCD prediction and the cor- 
rected data distribution, while <7y is the covariance matrix taking into account the sta- 
tistical error of the data, the finite Monte Carlo statistics, and the correlations between 
the bins of the distributions. The data and Monte Carlo event samples were each divided 
into 90 subsamples^ to compute a^, using the standard formula: 



1 

where Z)j is the value of the corrected distribution in bin i of the nth subsample, is the 
value in that bin for the whole sample and N is the number of subsamples. We included 
bins in the fit only if the total correction factors were in the range 0.9 < C\ ot < 1.1, except 
for the four-jet rate, where 0.85 < Cj tot < 1.15 was used.Q This bin choice is motivated 
by our expectation that the Monte Carlo estimate is less reliable if hadronization and/or 
detector distortions are large. The total number of bins used in the fit was 82, leaving 79 
degrees of freedom since there are three fitted parameters (77, x and y, see Section 3.5). 

3.5 Fit results 

The central values of the fit results for the parameter 77 and colour factor ratios x and y 
are listed in the upper section of Table 0. The renormalization scale was fixed at x^ — 1. 
The errors and correlations (p) in the table are purely statistical. 

We observe a strong correlation between x and y. We tested other variables (for 
instance, the three-jet rate) in an attempt to decrease these correlations, but did not 
observe a significant improvement. The correlations were taken into account when the fit 
results were converted to the standard QCD parameters a s (Mz) and colour factors using 
the definitions in Eqs. (|10|) and (12) and T R — 1/2. These latter results are given in the 
lower section of Table |2|. 



4 Systematic uncertainties 

The systematic uncertainties were evaluated by considering the following effects: 

• the measurement process and accuracy of the Monte Carlo detector simulation, 

• dependence on the model of hadronization, 

3 The number of subsamples is a result of optimization for the given data and number of bins. 
4 There are not enough bins in the smaller range to obtain a reliable fit. 
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rj = a s C F /(2n) 


0.0256±0.0003 


x = C A /C F 


2.25±0.08 


V = TnlC F 


0.37±0.04 


Prj x 


-0.33 


Pvv 


-0.11 


Px y 


0.90 


X 2 /d.o.f. 


98.5/79 


a s 


0.120±0.011 


C A 


3.02±0.25 


c F 


1.34±0.13 



Table 2: Results of the default fit. The errors and correlations are statistical only. 

• theoretical uncertainties arising from the choice of the renormalization scale, 

• variation of the matching scheme in the case of the differential two-jet distribution, 

• variation of the fit range, 

• the background from five parton events, 

• variation of the parameter y cnt . 

In the systematic checks we always used the same covariance matrix that was optimized 
for the standard analysis. 

4.1 Detector effects 

The uncertainty related to the detector simulation and the measurement process was 
estimated by repeating the analysis using either tracks only or clusters only for the data 
and the detector level Monte Carlo samples. For the tracks and clusters we applied 
the same selection procedure as described in Section 3.2, except that all cluster related 
criteria were left out in the former case and all track related criteria were left out in the 
latter. Furthermore, the entire analysis was repeated using | cos thrust! < 0.7 (rather than 
| cos thrust | < 0.9), or N c h > 8 (rather than N c h > 5). The results obtained for these 
analyses are given in Table [| 





rj x y x 2 /d-o.f. 


tracks only 
clusters only 
1 cos6» thrust | < 0.7 
N ch > 8 


0.0257±0.0003 2.24±0.08 0.37±0.04 127.1/79 
0.0256±0.0003 2.25±0.08 0.38±0.04 124.2/79 
0.0256±0.0003 2.26±0.08 0.38±0.04 117.9/79 
0.0256±0.0003 2.25±0.08 0.37±0.04 98.3/79 


standard analysis 


0.0256±0.0003 2.25±0.08 0.37±0.04 98.5/79 



Table 3: Fit results obtained by varying systematic conditions related to detector effects. 
The errors are statistical only. 
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4.2 Hadronization models 



Several combinations of different Monte Carlo models were considered to estimate the 
uncertainty associated with the hadronization correction. Since we used two different 
Monte Carlo samples to determine the correction factors (ME and PS), we examined the 
effects induced by variations of these models separately. 

The angular correlations were corrected using the ME version of JETSET with the 
Lund string fragmentation scheme, as stated in Section 3.3.1. In the systematic studies 
we used the ME version of the HERWIG generator |!9|, which produces four-parton 



final states according to the leading order QCD matrix elements, and then a parton 
shower is started from each of the four partons followed by hadronization based on the 
cluster fragmentation.^ We also varied the principal hadronization parameters of the 
JETSET ME Monte Carlo program, a q and a. As the hadronization parameter b is 
strongly correlated with a we kept it fixed. We increased a or decreased a g , alternately, by 
10 %. These changes correspond to about one standard deviation in the parameter values 
allowed for the JETSET ME generator at Z peak, as found by the LEP experiments |4(J . 
The fits were redone using data corrected by the factors calculated from these Monte 
Carlo samples at the generator level (reevaluating C] iad only) and using the standard Cf et 
correction factors for detector effects. 

The effect of the hadronization process on the four-jet rates and differential two-jet 
rates was studied by using the PS version of HERWIG instead of the PS version of JET- 
SET and by varying the Lund hadronization parameters a q and a as for the angular 
correlations. The parameter Q Q defining the end of the parton shower cascade was also 
varied from its standard value of 1.9 GeV to 2.4 GeV. The theoretical predictions assume 
zero quark masses. To assess the effect of this assumption, we determined the hadroniza- 
tion correction, C had , using JETSET Monte Carlo events with only light primary quarks 
(u, d, s, c) at the parton level but with all available flavours (u, d, s, c, b) at the hadron 
level. The results provided by these variations are presented in Table f|. 



4.3 Renormalization scale ambiguity 

The arbitrariness in the choice of the renormalization scale \i leads to an ambiguity in the 
theoretical prediction. We perform our default analysis at the scale factor x^ = \ij E cm = 1. 

To estimate the systematic error due to the scale ambiguity, we varied the renormal- 
ization scale within the range 0.5 < x^ < 2. Varying x^ probes the effect of missing 
higher orders, and these could be very different for the NLO and NLLA calculations. 
We performed two variations separately: varying x^ for the jet-related variables while 
keeping it fixed at the default value for the four-jet angular correlations, or vice versa. 
Table | presents the results of this systematic check. As expected, the normalised angular 
correlations were not sensitive to the scale choice, but the jet related variables showed 
significant sensitivity. In Fig. [3] we show the dependence of the x 2 °f the fit on x^, for 
the case of jet-related variables with x^ — 1 kept fixed for the angular correlations. The 
minimum of the x 2 is seen to lie close to i M = 1. 

5 The HERWIG-ME has not been tuned to OPAL data. We used CLSMR(1)=0.4, PSPLT(2)=0.33, 
DECWT=0.7, and parameters not listed here were left at their default values. 
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HW-ME+HW-PS 
JT-ME(a + la)+JT-PS* 
JT-ME(o-g - la)+JT-PS* 
JT-ME*+JT-PS(a + la) 
JT-ME*+JT-PS(a g - la) 

JT-ME*+JT-PS(Q + l(r) 
JT-ME*+JT-PS(wcfec) 



V 

0.0245±0 
0.0258±0 
0.0257±0 
0.0253±0 
0.0253±0 
0.0252±0 
0.0260±0 



.0003 
.0003 
.0003 
.0003 
.0003 
.0003 
.0003 



x 

2.43±0.18 
2.15±0.08 
2.16±0.08 
2.28±0.08 
2.28±0.08 
2.29±0.08 
2.19±0.07 



y 

0.44±0.08 
0.33±0.04 
0.33±0.03 
0.37±0.04 
0.37±0.04 
0.37±0.04 
0.33±0.03 



xVd-o.f. 
274.5/79 
94.5/79 
156.5/79 
106.6/79 
106.1/79 
101.1/79 
173.1/79 



standard analysis 



0.0256±0.0003 2.25±0.08 0.37±0.04 98.5/79 



Table 4: Fit results obtained by varying the Monte Carlo models used for the hadroniza- 
tion corrections. The abbreviations are as follows: * denotes the Monte Carlo model used 
in the standard analysis, the model before '+' is used to define hadronization corrections 
for the angular correlations while the model after the '+' is used to define hadronization 
corrections for the jet rates. JT-ME: JETSET 7.4 matrix element generator with string 
fragmentation, JT-PS: JETSET 7.4 parton shower generator with string fragmentation, 
HW-ME: HERWIG 6.1 matrix element generator with cluster fragmentation, HW-PS: 
HERWIG 5.9 parton shower generator with cluster fragmentation. 



4.4 Matching ambiguity 

For the differential two-jet rates, we performed the fit using the i?-matching scheme to 
combine the NLL and NLO approximations rather than the In i?-matching. The result of 
this variation is listed in Table [5[ 



4.5 Fit range 

The uncertainty related to the number of bins selected for the fit was evaluated by repeat- 
ing the fit using those bins for which the values of the correction factors were in the range 
of 0.8 < Ci tot < 1.2 instead of the default range 0.9 < C; tot < 1.1 or 0.85 < C; tot < 1.15. 
The results are presented in Table [5[ The variation of the fit range for the jet-related 
variables does not influence the value of the x 2 significantly. The fit-range variation for 
the angular variables cosxbz and cos is the main source of the much increased value 
of the x 2 - 



4.6 Five-parton background 

Five-parton events which are reconstructed as four-jet events at the parton level are 
taken into account by the NLO theory. We examined the effect of those five-parton 
events which yield four jets at the detector level but five at the parton level. To reduce 
this five-parton background, we applied a cut on y 45 , the value of the resolution variable 
which separates the four- and five-jet configurations of any given event. We optimized the 
value of this cut by studying five-parton events generated with the PENTAJET Monte 
Carlo program ||41|| . This program generates a five-parton configuration using the same 
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value of the intrinsic y$ as the JETSET-ME program, then hadronizes the partons using 



JETSET string fragmentation with the ERT-MC-3 parameters of Ref. |34]. The events 



were passed through the OPAL detector simulation program. Of 20 000 five-parton events, 
2 517 were found to be four-jet events at detector level, but five-jet events at parton level, 
for y cnt = 0.008. The number of these events reduces to 376 when we keep events with 
?/45 < 0.003, eliminating 85 % of the five-parton background, while leaving 60 % of the 
data. 

The results of the fit with this additional requirement are presented in Table [|. 



4.7 Dependence on the parameter y cut 

The value of the jet resolution parameter y cut , used in selecting four-jet events to define 
the angular correlations, can be chosen arbitrarily, as long as it is safely larger than the 
intrinsic yo of the JETSET-ME program. In order to have a high statistics sample, our 
default choice was the smallest possible value, which allows us to examine the stability of 
the results when y cut is increased. The four-jet event sample decreases very rapidly with 
increasing y cut , which constrains the region in y CVLt where a reliable systematic check can 
be performed. We repeated the entire analysis using y cnt = 0.01. The results of this fit 
are presented in Table [5| 





1] x y xVd.o.f. 


x M = 1.0 for angular corr. and 
= 0.5 for i? 4 and D 2 
= 2.0 for i? 4 and D 2 


0.0235±0.0003 2.40±0.08 0.28±0.03 118.0/79 
0.0271±0.0003 2.14±0.08 0.42±0.04 103.7/79 


x^ = 1.0 for i?4 and D 2 and 
x M = 0.5 for angular corr. 
x^ = 2.0 for angular corr. 


0.0256±0.0003 2.25±0.08 0.37±0.03 98.2/79 
0.0257±0.0003 2.24±0.08 0.37±0.04 99.0/79 


R-matching for D 2 


0.0281±0.0005 2.06±0.09 0.40±0.03 103.4/79 


0.8 < d tot < 1.2 


0.0255±0.0002 2.28±0.06 0.39±0.03 339.6/96 


5-parton background 


0.0257±0.0003 2.18±0.08 0.34±0.04 112.5/79 


ycut = o.oi 


0.0257±0.0003 2.27±0.21 0.38±0.09 127.1/79 


standard analysis 


0.0256±0.0003 2.25±0.08 0.37±0.04 98.5/79 



Table 5: Fit results obtained by varying the renormalization scale, the matching scheme 
for D 2 , the fit range, the 5-parton background and the y cut parameter. 



4.8 Additional systematic checks 

We tested the influence of the two- and three-parton events that give rise to four-jet events 
at the detector level by repeating the analysis requiring E m [ n = 5 GeV on the minimum 
energy of the jets instead of the default value E m i n = 3 GeV. A higher value of the jet 
energy reduces the number of two- and three-parton events which satisfy the four-jet 
selection criteria in the detector. We generated 1 million events with two, three and four 
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partons in the initial state using the ME version of JETSET with renormalization scale 
set to x 2 ^ = 0.002, which provides the correct proportion of two, three and four jet events 
for y cut = 0.008. We then counted the proportion of the two- and three-parton events 
which passed the four-jet selection at detector level. With E > E min = 5 GeV, there are 
no background events from the two- and three-parton events. With E > E m i n = 3 GeV, 
a background of 1 % appears (519 two- or three-parton events in the sample of 53 302 
four-jet events selected at the detector level). We found negligible effects on the final 
results of the fit, therefore, we did not include the effect of this check in the systematic 
errors. 

As stated in Section 2, the fits to individual angular correlations give rather ambiguous 
results. Therefore, we choose to use all four angular correlation variables in our standard 
analysis. As a consistency check, we performed the fit using only five variables instead 
of six, removing one of the six variables at a time. The results of these fits are presented 
in Table |6|. The errors are statistical only. We found that leaving out any one of the 
observables yields consistent results for the measured parameters within the statistical 
errors of the measurements. Therefore, we did not include the effect of this check in the 
systematic errors. 

The number of light quark flavours in a fixed order calculation is usually set to Nf = 5 
at LEP energies because the hard scattering scale is assumed to be the e + e~ center of 
mass energy. However, one may argue that in the four-quark channel the g —>■ bb splitting 
is suppressed compared to the g — ► qq (q=u, d, s, c) splittings, because the scale of 
that process is much lower than the scale at the primary vertex. To study the effect of 
suppressed g — > bb splitting on our results, we also used Nf = 4 light flavours at the 
secondary vertex in the fixed order prediction for the angular correlations. The results of 
this check are also shown in Table ^. The colour factor ratios increase slightly, remaining 
well within the statistical errors. The quality of the fit is almost unchanged. We did not 
include the effect of this check in the systematic errors, because it is not customary to do 
so. 





r] x y x 2 /d.o.f. 


without 




Bengtsson-Zerwas angle 


0.0257±0.0003 2.22±0.17 0.36±0.08 63.7/64 


Nachtmann-Reiter angle 


0.0256±0.0003 2.39±0.20 0.44±0.09 56.0/62 


Korner- S chierholt z- Willro dt angle 


0.0256±0.0003 2.29±0.21 0.39±0.09 63.1/63 


angle of the two softest jets 


0.0256±0.0003 2.31±0.17 0.40±0.08 57.9/65 


four-jet rates 


0.0243±0.0007 2.18±0.09 0.39±0.04 75.8/73 


differential two-jet rates 


0.0251±0.0010 2.31±0.12 0.36±0.04 52.1/65 


N f = A 


0.0256±0.0003 2.28±0.08 0.39±0.04 95.9/79 


standard analysis 


0.0256±0.0003 2.25±0.08 0.37±0.04 98.5/79 



Table 6: Fit results using five variables only or using Np = 4. 
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4.9 Total systematic error 

To combine the effects of the systematic variations we employed a Bayesian method that 
was developed in a similar analysis by the ALEPH collaboration M, where the explicit 
formulae we used can be found. The basic idea is that the information obtained from a 
fit is accepted or rejected based on the quality of the fit, i.e. the magnitude of the x 2 . 
From a Bayesian point of view we presume that all models are equally well suited to the 
analysis. However, a large x 2 /d.o.f. indicates that the probability of a model is low and, 
therefore, it should contribute to the systematic errors with a small weight. Accordingly, 
the systematic error corresponds to an increase of the x 2 by 1. 

We could not use the Bayesian method in two cases. In the variation of the fit range the 
number of bins changed, while in checking the renormalization scale uncertainty for the 
angular correlations we did not find a pronounced minimum of the x 2 curve. Furthermore 
we did not use the Bayesian method to estimate the systematic uncertainty from the five- 
parton background because the event statistics for this check were modified significantly, 
by about 40%, relative to the standard analysis. In these cases we estimated the systematic 
error by the deviation of the results from the standard values. 

The errors obtained from the different systematic checks are summarized in Table 0. 
The theoretical uncertainty refers to the larger of the errors from the scale uncertainty 
of the jet-related variables, the matching ambiguity (both of those estimate the effect of 
unknown higher order perturbative contributions), and the scale uncertainty of the four- 
jet angular variables. The different sources were added in quadrature to define the total 
systematic error .0 

Taking into account the systematic correlations presented in Table [7|, we converted 
these values to the errors of the standard QCD parameters a s (Mz) and colour factors 
using the definitions in Eq. (|T0|) and (|T2|) and Tr = 1/2. The result of the conversion is 
presented in Table [7]. The errors are dominated by the theoretical uncertainty. 





Ai] Ax Ay p vx p m p xy 


Detector effects 


0.0000 0.004 0.002 -0.87 -0.89 0.95 


Hadronization 


0.0003 0.056 0.021 -0.86 -0.71 0.92 


Theoretical uncertainty 


0.0013 0.096 0.038 -1.0 1.0 -1.0 


Fit range 


0.0001 0.033 0.013 -1.0 -1.0 1.0 


5-parton background 


0.0001 0.070 0.039 -1.0 -1.0 1.0 


Vcut 


0.0000 0.005 0.002 1.0 1.0 1.0 


Total systematic uncertainties 


0.0013 0.136 0.060 -0.83 0.06 0.38 




Aa s AC A AC f 


Total systematic uncertainties 


0.020 0.49 0.22 



Table 7: Contributions to the total systematic error. 



6 If the Bayesian approach were not adopted and the systematic errors were estimated by the largest 
deviation from the standard values, then we obtain the following systematic errors: Arj = 0.0027, Ax — 
0.27, Ay = 0.08. For the basic QCD parameters these errors are: Aa s = 0.030, AC A = 0.72, ACf = 
0.31. 
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5 Results and Discussion 



The procedure outlined above yields 
Ca/C f = 2.25 ± 0.08(stat.) ± 0.14(syst.) , T R /C F = 0.37 ± 0.04(stat.) ± 0.06(syst. 

Repeating the entire analysis with a s (M z ) fixed to its world average value 0.1184 
leads to similar results: 

C A /C F = 2.29 ± 0.06(stat.) ± 0.14(syst.) , T R /C F = 0.38 ± 0.03(stat.) ± 0.06(syst.) . 

To compare the results to previous measurements, we show in Fig. f| the two dimensional 
68% and 95% confidence level contour plots of the colour factor ratios based on total 
uncertainties obtained by our analysis. 

The DELPHI Collaboration ]I2| has performed a least-squares fit of leading order 
predictions to the two-dimensional distribution in the variables cosOnr and cos 034 in 
order to measure of the colour factor ratios. The DELPHI Collaboration [4"2"| also per- 



formed this analysis with b-tagging to separate quark jets from gluon jets which increases 
sensitivity. The OPAL Collaboration fitted the LO prediction to a three-dimensional dis- 
tribution of the cosGnr , cosa 34 and cos^bz variables flBfl . In a second measurement 
the OPAL Collaboration used event shape distributions of LEP1 data for fitting a s and 



one of the colour factors ||43|| . Furthermore, another recent measurement of the QCD 
colour factors and a s using event shape distributions at several energy points and power 
corrections was carried out ||44|| . The ALEPH Collaboration || was the first to present 



a simultaneous measurement of the strong coupling and the colour factors. The distribu- 
tion D2 and all four angular distributions defined in Section 2 have been employed. The 
results of the ALEPH experiment were determined using leading order predictions for the 
four-jet observables. The systematic uncertainties include contributions from variation 
of the renormalization scale, the matching ambiguity, from hadronization and detector 
effects, and finally from the estimation of mass effects. Our results are in agreement with 
these previous results. The main new feature of our analysis is that we obtained our re- 
sults using NLO predictions for four-jet angular correlations. Furthermore, our systematic 
checks contain additional items arising from the estimate of the effects of the five-parton 
background and the variation of the y cut parameter. We also checked that assuming four 
rather than five light flavours at the secondary vertex does not influence our results. This 
effect was not examined in any of the previous studies. 

With the inclusion of the higher order corrections, the colour factor ratios became 
better constrained, especially T R /C F . We found an increase of about 15% on T R /C F 
using the NLO theoretical predictions instead of the LO one, while Ca/C f was only 
slightly affected (about 3%). 

The diamond symbol in Fig. (| indicates the result one would expect if a leading order 
perturbative theory with light gluinos was the correct theory. This point falls well outside 
our 95% C.L. contour, indicating that the data do not favour the existence of light gluinos. 
A complete analysis to test for the existence of gluinos should use hadronization models 
and NLO predictions that include the effect of the gluinos 0. The latter is possible using 
the DEBRECEN partonic Monte Carlo program we used for generating the NLO results (see 
the Appendix), but presently none of the existing hadronization models include gluinos. 
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Finally, we converted our measured parameters to the standard QCD parameters 
(using T R = 1/2), which leads us to our main results: 



a s = 0.120 ± 0.011(stat.) ± 0.020(syst.) 
C A = 3.02 ± 0.25(stat.) ± 0.49(syst.) 
C F = 1.34 ± 0.13(stat.) ± 0.22(syst.) 



Other techniques |22| provide precise determinations of a s . In this study the measure- 



ment of a s serves as a crosscheck and a consistent value was found. The larger uncertainty 
we obtain for a s compared to other techniques can be traced to the functional dependence 
on the other two parameters. 

6 Summary 

In this paper, a test of perturbative QCD at LEP has been presented. With only the 
assumption of non-abelian gauge symmetry and standard hadronization models, we mea- 
sured the eigenvalues of quadratic Casimir operators — the colour factors — of the un- 
derlying gauge group, together with the coupling of QCD. The measurement was based 
on the comparison of next-to-leading order perturbative predictions for four-jet angular 
variables, and the resummation improved next-to-leading order perturbative predictions 
for the multi-jet rates, to OPAL data taken between 1991 and 1995 at the Z peak. The 
results for the colour factors follows: 

C A = 3.02 ± 0.25(stat.) ± 0.49(syst.) , C F = 1.34 ± 0.13(stat.) ± 0.22(syst.) , 

from a global fit to the angular correlations, four-jet rates, and differential two-jet rates. 
These values are in agreement with SU(3) values of Ca = 3 and Cf = 4/3 as well as with 
measured values obtained previously at LEP. The corresponding strong coupling value, 

a s = 0.120 ± 0.011(stat.) ± 0.020(syst.) , 

is in agreement with the world average. 



A Appendix: Theoretical Predictions 

The NLO perturbative expansion of the differential two-jet rates has the following expres- 
sion: 



D 2 {y 



1 da 



23 I 



otot dy ; 



2:-! 



v(v)A' D2 (y 23 ) + V (fi) 2 \B' D2 (y 23 ) + f3 \n(x^)A' D2 (y 23 ] 



where y 23 is the y cnt value which separates the two- and three-jet configurations of an 
event. The corresponding formula in the case of the four-jet rates takes a similar form: 



cut J 



Ctot 



cut , 



(9) 
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In Eqs. (H) and (Q) cr to t is the total hadronic cross section at 0(r/) accuracy and 



oJjACf r,(Mz) /, A ... , \n(w(^,M z )) \ 

riip) = ~^r~ = iP3 V ~ s ,( z) u^,M z) j • (10) 

In the second equality in Eq. (|T0|), we used the two-loop expression for the running 
coupling with 

w(ji,M z ) = l-Pa7,(jJL)]n(^j , (11) 

and \i denotes the renormalization scale (x M = fi/y/s, where y/s is the total cm. energy). 
We define the colour factor ratios as 

T r 

X = -FT, y = TT, Vf = N f y, (12) 

where Nf is the number of active flavours. Then the coefficients in the expansion of the 
j3 function can be written as 

11 4 . 17 2 10 

Po = yX --?//, fi l = —x -2y } - —xy f . (13) 

The functions A' D2 , Br 4 and B' D2 , Cr 4 in Eqs. @ and @ are the perturbatively calculable 
coefficient functions in the Born approximation and the radiative correction, respectively, 
which are independent of the renormalization scale. The A' D2 function is independent of 
x and y. The B functions are linear forms and the C4 function is a quadratic form of 
these colour charge ratios 



where 



with 



B = B { ° ] + BP x + B^y, O = D 2 , or R 4 (14) 
C R4 = + C« x + CW y + C« z + C« x 2 + Cg> * y + C« , (15) 



iV c C| 



(16) 



C c = ^ Tr(T a T 6 T tc )Tr(T tc T b T a ) (17) 



a,ft,c=l 



being a cubic Casimir invariant of the gauge group.[] We calculated these coefficient 
functions of the differential two-jet rates in the —6 < 01(3/23) < 2/23 range and those of 
the four-jet rates in the — 3 < log 10 (y cu t) < y cnt range using the partonic Monte Carlo 
program DEBRECEN [ 46| . All theoretical predictions were obtained for Nf = 5 light-quark 
flavours and with the normalization Tr = |. The latter choice is arbitrary and influences 
the values of the coefficient functions in such a way that the physical prediction remains 
unchanged. 

7 The results of the Monte Carlo integration for ci 4 ' are compatible with zero [pi, therefore, z is not 
included as a free parameter in our fits. 
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The resummation formulas for various event shapes are given in Ref. T? at the NLL 
accuracy for the cumulative cross section, 

fy 1 da 

R(y)= / ch/ 23 -f^, (18) 
Jo a dy 23 

in the following form: 

R NLL (L) = (1 + drj + c 2V 2 )e L9l ^ +92 ^ , (19) 

where L = — ln(y 23 ). The functions Lgi(rjL) and g 2 (rjL) represent the sums of the leading 
and next-to-leading logarithms, respectively, to all orders in the strong coupling. The 
colour decomposition of the g± and g 2 functions is given in Ref. The coefficient 



C\ = — 5/2 + 7r 2 /6 — 6 In 2 is known exactly ||9| and C 2 can be determined from numerical 
integration of the 0(a%) matrix elements (see below). 

At NLO (i.e. O(aig) accuracy), the cumulative cross section for the differential two-jet 
rates has the form: 

R 0(aj) (L) = 1 + A(L) 7] + B(L)r] 2 . (20) 

In order to have the best possible theoretical description, we fit the experimental data 
to the matched NLL and NLO results, which is expected to describe the data in a wider 
kinematical range than any of the two separately. In the case of D 2 , our default match- 
ing procedure is the In i?-matching, which combines the NLL and NLO approximations 
according to the following formula: 



+A(L) rj + [ B (L) + x B X {L) + y B y (L) - -A(L) 2 ) r, 2 



In R(L) = E(L, rj) - (G n L + G 12 L 2 ) V - (G 22 L 2 + G 23 L 3 ) V 2 (21) 

1 

2" 

where E(L, rj) = Lgi(rjL) + g 2 (rjL) and the following two terms represent the two lowest 
order terms in the ^-expansion of InR. In the actual analysis, we used the differential 
distribution da/dL that is obtained from Eq. (|2T| ) as 



1 da ^. 
-— = exp(lni?) 
a dL 



dS(L, 77) _ + 2 ^ ^ _ ^2 G 22 L + 3 G 23 L 2 ) V 2 (22) 



dL 

+A'(L) v + (B' (L) + x B' X (L) +yB' y (L) - A(L) A\L)) rf 

with A'(L) and B'(L) given in Eq. (@). 

Eq. (|22]) contains the cumulative coefficients A(L) and Bi(L), the latter ones implicitly 
in R, that cannot be directly calculated using a partonic Monte Carlo program. We know, 
however, that 

A(L) = C 1 + y £ R lm L m + D x (L) , (23) 

m=l 

and 

Bi{L) = Cf + J2 R$l L m + Df{L) , i = 0,x,y, (24) 

m=l 
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where the colour decomposition of the R^m logarithmic coefficients can be found in 
Ref. H, except for R 2 l- This latter coefficient as well as the C% constants and the 
Di(L), Df{L) functions can be determined by differentiating Eqs. ( p3|) and (^4j) with 
respect to L and fitting to the coefficient functions as obtained from the partonic Monte 
Carlo program DEBRECEN. Following the fitting procedure of Refs. [[IT], [48] we obtained 
ff§} = -16.531 ± 2.97(stat.), R$ = -4.293 ± 1.03(stat.), R$ = 4.54 ± 1.05(stat.), 
C 2 ° = -25.56 ± 4.04(stat.), Cf = -1.53 ± 1.12(stat.), Cf = -3.565 ± 1.17(stat.). 

In order to account for the dependence on the renormalization scale x^, one has to 
make the following substitutions in Eq. (|22|) : 

B(L) ^ B(L) + A(L) {3 \nx^ (25) 
E(X,r 1 )^E(X, V ) + X 2 g' 1 \nx IM , (26) 
G 22 ^G 22 + G 12 f3 \n Xfl . (27) 

The other widely used matching scheme is i?-matching, which will serve as a systematic 
check of our matching procedure. In this case the NLL and NLO predictions are combined 
according to the following formula |p],[l6[]: 

^R-match = ^NLL + U ( B - £ NLL ) + 7] 3 (C - C NLL )1 , (28) 



where _B NLL and C NLL are the coefficients in the ^-expansion of i? NLL . In the case of multi- 
jet rates, the NLL approximation does not exponentiate, therefore, the viable matching 
scheme for the four-jet rates is i?-matching. The resummation formulas for the four-jet 



rates can be found in Refs. 16. 18 



The general form of the NLO differential cross section for the four-jet angular corre- 
lations O4 is given by |TJj : 



--^r(O i )= V ( f i) 2 B 0i (0 4 )+r ] ( f i) 3 [CoAOa) + 5o 4 (0 4 ) A) H^)} , (29) 

Cq O.U4 L J 

where ctq denotes the Born cross section for the process e + e~ — > qq. In Eq. (|29|) the Born 



and correction functions B 0i and Co 4 have an analogous decomposition to Eqs. (fTJ) 
and fllpp. To obtain the distributions normalised to unity in perturbation theory, we 
calculated the coefficient functions b\°^ and C\°^ in the linear and quadratic forms 
using the partonic Monte Carlo program DEBRECEN, and normalised those with the total 
0{rf) four-jet cross section of the same y cut . 
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Figure 1: Distributions of the four-jet rates and the differential two-jet rates. OPAL data 
corrected to the parton level are denoted by points. The statistical uncertainty in the 
bins is smaller than the size of the circles. The small sections below each plot indicate the 
correction factors C\ ot for experimental and hadronization effects with their statistical 
error (shaded area). The horizontal lines with arrows indicate the selected fit region. 
The QCD predictions shown by the dotted and dashed lines employ the world average 
a s =0.119 and the standard values Ca = 3 and Cp = 4/3. The solid lines represent the 
result of the simultaneous fit of all six variables. 
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Figure 2: Distributions of the four-jet angular correlations. OPAL data corrected to the 
parton level are denoted by points. The statistical uncertainty in the bins is smaller 
than the size of the circles. The small sections below each plot indicate the correction 
factors C\ ot for experimental and hadronization effects with their statistical error (shaded 
area). The horizontal lines with arrows indicate the selected fit region. For the normalised 
angular correlations, the distributions at leading and next-to-leading order are almost the 
same 



14, 19 1, therefore, only the NLO results are shown. 
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Figure 3: Dependence of the fitted results for the colour factor ratios and r\ = a s Cjr/(27r) 
on the renormalization scale factor x^. The shaded area represents the range allowed by 
the statistical errors. The horizontal lines show the SU(3) values and the world average 
for a s (Mz) translated to rj using Eq. (0). The bottom fi sure shows the variation of the 
X 2 of the fit with x M varied for jet-related variables. 
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Figure 4: Two dimensional plot of the ratios of the colour factors compared with previous 



results of the LEP experiments P, |T2j, |13| . The legend in the top right corner describes 
this analysis. The contours are based on total uncertainties. 
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